home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / CUJ9204.ARJ / 1004024A < prev    next >
Text File  |  1992-06-02  |  4KB  |  129 lines
  1. /* Listing 1 */
  2. typedef struct {
  3.     double X1, X2, X3, X4;
  4. }      ARG_D_4;            /* vector 4 */
  5. #include "float.h"
  6. ARG_D_4 sin_4(xx1, xx2, xx3, xx4)
  7.     double xx1, xx2, xx3, xx4;
  8. /*  use where cos of same argument not needed
  9. ** 16 digits precision, compare to 15 digits in "dtrig.h"
  10. ** T C Prince */
  11. {
  12.     union dblfmt {
  13.     double dbl;
  14.     int idbl;
  15.     struct dfmt {        /* IEEE p754 */
  16.         unsigned int sign:1;
  17.         unsigned int ex:11;
  18.     }    fmt;
  19.     }      xi1;
  20.     double xr, x2, x4, x8;
  21. #ifdef __STDC__
  22.     long double pi = 3.1415926535897932385, pml;
  23. #include <limits.h>
  24. #else
  25.     register double pi = 3.1415926535897932385, pml;
  26. #define LONG_MIN 0x80000000
  27. #endif
  28.     union dblfmt pm, round;
  29.     ARG_D_4 res;
  30. #define BIAS DBL_MAX_EXP
  31. #include <errno.h>
  32. #ifndef errno
  33.     extern int errno;
  34. #endif
  35. #define ODD(i) ((i)&1)
  36. /* use identity sin(x + n pi) = (-1)^n sin(x)
  37. **  to reduce range to -pi/2 < x < pi/2
  38. **    pml=rint(xi1/pi) */
  39. #if FLT_ROUNDS != 1
  40.     #error "rounding mode not nearest; adjust code"
  41. #endif
  42. #if FLT_RADIX !=2 && FLT_RADIX != 10
  43.     #error "code not optimum for accuracy in this RADIX"
  44. #endif
  45. #if DBL_DIG > 16
  46.     #error "more terms needed for full accuracy"
  47. #endif
  48. /* shortcut test of sign, not portable to VAX */
  49.     round.dbl = 1 / LDBL_EPSILON;
  50.     xi1.dbl = xx1;
  51.     round.idbl |= xi1.idbl & LONG_MIN;
  52.     pml = xx1 / pi + round.dbl;
  53. /* sign reversal may reduce register usage */
  54.     xr = pi * (pml -= round.dbl) - xx1;
  55. /* shortcut test for fabs(pml) > INT_MAX */
  56.     pm.dbl = pml;
  57.     if (pm.fmt.ex > BIAS + 31)
  58.     errno = ERANGE;
  59. /* don't wait to calculate xr**2 until sign is fixed;
  60. ** another sign reversal is due if pm.dbl is odd */
  61.     x2 = xr * xr;
  62. /* first sign reversal compensated in coefficient signs;
  63. ** conditional sign fixed by testing odd/even
  64. ** first two results are obtained by straight Horner
  65. ** polynomial evaluation */
  66.     res.X1 = (-.9999999999999999 + x2 * (.1666666666666607
  67.     + x2 * (-.833333333328281e-2 + x2 * (.19841269824875e-3
  68.     + x2 * (-.2755731661057e-5 + x2 * (.25051882036e-7
  69.      + x2 * (-.160481709e-9 + x2 * .7374418e-12)))))))
  70.     * (ODD((int) pm.dbl) ? -xr : xr);
  71. /* sin(xi2) */
  72.     round.dbl = 1 / LDBL_EPSILON;
  73.     xi1.dbl = xx2;
  74.     round.idbl |= xi1.idbl & LONG_MIN;
  75.     pml = xx2 / pi + round.dbl;
  76.     xr = pi * (pml -= round.dbl) - xx2;
  77.     pm.dbl = pml;
  78.     if (pm.fmt.ex > BIAS + 31)
  79.     errno = ERANGE;
  80.     x2 = xr * xr;
  81.     res.X2 = (-.9999999999999999 + x2 * (.1666666666666607
  82.     + x2 * (-.833333333328281e-2 + x2 * (.19841269824875e-3
  83.     + x2 * (-.2755731661057e-5 + x2 * (.25051882036e-7
  84.      + x2 * (-.160481709e-9 + x2 * .7374418e-12)))))))
  85.     * (ODD((int) pm.dbl) ? -xr : xr);
  86. /* sin(xi3) */
  87.     round.dbl = 1 / LDBL_EPSILON;
  88.     xi1.dbl = xx3;
  89.     round.idbl |= xi1.idbl & LONG_MIN;
  90.     pml = xx3 / pi + round.dbl;
  91.     xr = pi * (pml -= round.dbl) - xx3;
  92.     pm.dbl = pml;
  93.     if (pm.fmt.ex > BIAS + 31)
  94.     errno = ERANGE;
  95.     x2 = xr * xr;
  96.     x4 = x2 * x2;
  97. /* split into 2 Horner polynomials to increase
  98. ** parallelism after 1st result finishes */
  99.     res.X3 = (-.9999999999999999 + x2 * (.1666666666666607
  100.                   + x2 * (-.833333333328281e-2
  101.                    + x2 * .19841269824875e-3))
  102.           + (-.2755731661057e-5 + x2 * (.25051882036e-7
  103.                     + x2 * (-.160481709e-9
  104.                + x2 * .7374418e-12))) * x4 * x4) *
  105.     (ODD((int) pm.dbl) ? -xr : xr);
  106. /* sin(xi4) */
  107.     round.dbl = 1 / LDBL_EPSILON;
  108.     xi1.dbl = xx4;
  109.     round.idbl |= xi1.idbl & LONG_MIN;
  110.     pml = xx4 / pi + round.dbl;
  111.     xr = pi * (pml -= round.dbl) - xi1.dbl;
  112. /* errno is set to ERANGE if any of the arguments are too
  113. ** large for reasonable range reduction */
  114.     pm.dbl = pml;
  115.     if (pm.fmt.ex > BIAS + 31)
  116.     errno = ERANGE;
  117.     x2 = xr * xr;
  118.     x4 = x2 * x2;
  119.     x8 = x4 * x4;
  120. /* multiply by 1 is K&R way to enforce parentheses */
  121.     res.X4 = ((-.9999999999999999 + x2 * (.1666666666666607
  122.                   + x2 * (-.833333333328281e-2
  123.               + x2 * .19841269824875e-3))) * 1
  124.           + (-.2755731661057e-5 + x2 * (.25051882036e-7
  125.        + x2 * (-.160481709e-9 + x2 * .7374418e-12))) * x8)
  126.     * (ODD((int) pm.dbl) ? -xr : xr);
  127.     return res;
  128. }
  129.